This code simulates the Keen model as specified in Grasselli and Costa Lima (2012)
Authors: B. Costa Lima and M. Grasselli
Grasselli, M. R., and B. Costa Lima. 2012. “An Analysis of the Keen Model for Credit Expansion, Asset Price Bubbles and Financial Fragility.” Mathematics and Financial Economics 6 (3): 191–210. https://doi.org/10.1007/s11579-012-0071-8.
"C:\Users\ibuckley\OneDrive - CSTO\Code\Grasselli\2012 Keen\"
# using Pkg;
# Pkg.add("ParameterizedFunctions")
# Pkg.add("DifferentialEquations")
# Pkg.add("PlotlyBase")
using DifferentialEquations # for ODEProblem
using DataFrames
using Plots
plotly() # Choose Plotly as the backend to Plots
┌ Info: Precompiling DataFrames [a93c6f00-e57d-5684-b7b6-d8193f3e46c0] └ @ Base loading.jl:1260
Plots.PlotlyBackend()
#using Plotly
#Plotly.signin("IanBuckley","Z4HJZpdeW4wss1tnizq8")
where the rate of new investment is a nonlinear increasing function $\kappa$ of the net profit share $π = (1 − ω − rd)$ and $δ$ is a constant depreciation rate as before.
The new dynamic variable in this model is the amount of debt, which changes based on the difference between new investment and net profits.
ν = 3 # capital-to-output ratio
α = 0.025 # productivity growth rate
β = 0.02 # population growth rate
δ = 0.01; # depreciation rate
ϕ0 = 0.04/(1-0.04^2) # Phillips curve parameter, Phi(0.96)=0
ϕ1 = 0.04^3/(1-0.04^2); # Phillips curve parameter, Phi(0)=-0.04
r = 0.03 # real interest rate κ
κ0 = -0.0065 # investment function parameter
κ1 = exp(-5) # investment function parameter
κ2 = 20; # investment function parameter
Φ(λ) = ϕ1./(1-λ).^2- ϕ0 # Phillips curve, equation (25)
Φ′(λ) = 2*ϕ1./(1-λ).^3 # derivative of Phillips curve
Φ_inv(λ) = 1 - (ϕ1/(λ+ϕ0)).^(1/2); # inverse of Phillips curve
κ(x) = κ0+κ1.*exp(κ2.*x) # Investment function, equation (64)
κ′(x) = κ1.*κ2.*exp(κ2.*x) # Derivative of investment function
κ_inv(x) = log((x-κ0)./κ1)./κ2; # Inverse of investment function
See Equation (42)
κ_eq = ν*(α+β+δ) # investment at interior equilibrium
π_eq = κ_inv(ν*(α+β+δ)) # profit at interior equilibrium, equation (41)
d_eq = (κ(π_eq)-π_eq)/(α+β) # debt ratio at interior equilibirum, equation (42)
ω_eq = 1-π_eq-r*d_eq # wage share at interior equilibrium, equation (42)
λ_eq = Φ_inv(α); # employment at interior equilibirum, equation (42)
# TO DO!
#d_sols(1) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,d) ; # solution to equation (40) given in equation (66)
#d_sols(2) = vpasolve(d*(r+delta-fun_kappa(1-r*d)/nu)+fun_kappa(1-r*d)-1,[20 40]) ; # solution to equation (40) given in equation (66)
#d_sols
$r\left[\frac{\kappa^{\prime}\left(\bar{\pi}_{1}\right)}{v}\left(\bar{\pi}_{1}-\kappa\left(\bar{\pi}_{1}\right)+v(\alpha+\beta)\right)-(\alpha+\beta)\right]>0$
J(𝜔,𝜆,𝑢)= [Φ(𝜆)-𝛼 𝜔.*Φ′(𝜆) 0;
-𝜆.*𝜅′(1-𝜔-r*z)/𝜈 (𝜅(1-𝜔-r*z)-𝜈*(𝛼+𝛽+𝛿))/𝜈 -r*𝜆.*𝜅′(1-𝜔-r*z)/𝜈;
((z-𝜈).*𝜅′(1-𝜔-r*z)+𝜈)/𝜈 0 (𝜈*(r+𝛿)-𝜅(1-𝜔-r*z)+r*(z-𝜈).*𝜅′(1-𝜔-r*z))/𝜈];
# TO DO!!!
#E1=eig(J(omega_eq,lambda_eq,d_eq)) # eigenvalues of the Jacobian at the interior equilibrium
#E2=eig(J(0,0,d_sols(1))) # eigenvalues of the Jacobian at the first d0 equilibrium
#E3=eig(J(0,0,d_sols(2))) # eigenvalues of the Jacobian at the second d0 equilibrium
cond59= κ′(π_eq)/ν*(π_eq-ν*δ)-(α+β) # expression in brackets in condition (59)
0.10573866662583055
DifferentialEquations¶function keen!(dz,z,p,t)
log_ω, tan_λ, π_n = z # state variables
α, β, δ, ν, r, (ϕ0,ϕ1), (κ0,κ1,κ2) = p # parameters & functions κ, Φ ??? ϕ0,ϕ1,κ0,κ1,κ2
# Recover λ & ω
λ = atan(tan_λ)/π + 0.5
ω = exp(log_ω)
# Functions
Φ(λ) = ϕ1./(1-λ).^2 - ϕ0 # Phillips curve, equation (25)
κ(x) = κ0 + κ1.*exp(κ2.*x) # Investment function, equation (64)
g_Y = κ(π_n)/ν - δ # Output function (31)
# Equation (47)
dz[1] = dlog_ω = Φ(λ) - α #d(log_ω)/dt
dz[2] = dtan_λ = (1 + tan_λ^2)*π*λ*(g_Y - α - β) #d(tan_λ)/dt (constant π!!)
dz[3] = dπ_n = -ω*dlog_ω - r*(κ(π_n) - π_n) + (1 - ω - π_n)*g_Y #d(π_n)/dt
end
keen! (generic function with 1 method)
convert¶function convert(state_variables,r)
n = length(state_variables)
if n==2
ω, λ = state_variables
elseif n==3
ω, λ, d = state_variables
else
ω, λ, d, p = state_variables
end
log_ω = log.(ω)
tan_λ = tan.((λ-0.5)*π) # π the constant!
if n>2
π_n = 1 - ω - r * d
if n==4
log_p = log.(p)
end
end
if n==2
log_ω, tan_λ
elseif n==3
log_ω, tan_λ, π_n
else
log_ω, tan_λ, π_n, log_p
end
end
convert (generic function with 1 method)
restore¶function restore(transformed,r)
n = length(transformed)
if n==2
log_ω, tan_λ = transformed
elseif n==3
log_ω, tan_λ, π_n = transformed
else
log_ω, tan_λ, π_n, log_p = transformed
end
ω = exp.(log_ω)
λ = atan.(tan_λ)/π + 0.5 # π the constant!
if n>2
d = (1 - ω - π_n)/r
if n==4
p = atan.(log_p)
end
end
if n==2
ω, λ
elseif n==3
ω, λ, d
else
ω, λ, d, p
end
end
restore (generic function with 1 method)
ODEProblem and solve¶# Sample paths of the Keen model
π_bar_n_K = κ_inv(ν * (α + β + δ));
d_bar_K = (κ(π_bar_n_K) - π_bar_n_K)/(α + β);
ω_bar_K = 1 - π_bar_n_K - r * d_bar_K;
λ_bar_K = Φ_inv(α);
ω0 = 0.75
λ0 = 0.75
d0 = 0.1
Y0 = 100
T = 300;
ODEProblem¶u0 = [-0.288, 1.0, 0.247]
tspan = (0.0,300.0)
p = [α, β, δ, ν, r, (ϕ0,ϕ1), (κ0,κ1,κ2)]
prob = ODEProblem(keen!,u0,tspan,p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 300.0) u0: [-0.288, 1.0, 0.247]
solve¶sol = solve(prob,reltol=1e-6);
DataFrame¶df = DataFrame(sol) # It's the wrong way around!!!!!!!
| x1 | x2 | x3 | x4 | x5 | x6 | x7 | x8 | |
|---|---|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | -0.288 | -0.293876 | -0.305295 | -0.320146 | -0.330904 | -0.33368 | -0.322919 | -0.312928 |
| 2 | 1.0 | 1.12342 | 1.46596 | 2.41279 | 4.98538 | 11.5144 | 35.2791 | 49.8736 |
| 3 | 0.247 | 0.24955 | 0.254226 | 0.259767 | 0.263179 | 0.262936 | 0.254095 | 0.246677 |
#state_variables = map(u -> restore(u,r),sol.u);
state_variables = [(restore(u,r),t) for (u,t) in zip(sol.u,sol.t)];
state_variables[2]
((0.7461811191854018, 0.7649302031586827, 0.15779923733286), 0.0748259153505046)
plot([u for (u,t) in state_variables])
xaxis!("ω - wage share of the economy")
yaxis!("λ - employment rate")
#zaxis!("d - debt")
[[u1,u2,u3] for ((u1,u2,u3),t) in state_variables]
1330-element Array{Array{Float64,1},1}:
[0.7497615922390413, 0.75, 0.10794692536529131]
[0.7461811191854018, 0.7649302031586827, 0.15779923733286]
[0.7398108038589221, 0.7945265193981843, 0.25112337052938993]
[0.7321778575222391, 0.8358964387697227, 0.37057465404153334]
[0.7263666965926021, 0.8727191304039716, 0.46742844682300577]
[0.7211036452946267, 0.9117273336598325, 0.5615264153944218]
[0.7182526371418824, 0.937183540173052, 0.618696007939749]
[0.7166849652319078, 0.9558377532286569, 0.6586181390481579]
[0.7162133455271785, 0.9690755840265715, 0.6859638726882704]
[0.716804590918236, 0.978543607812272, 0.7050051718226277]
[0.7187297827638306, 0.9853306659147664, 0.7183537329321847]
[0.7229203640146292, 0.9902674385001823, 0.7278498182572812]
[0.7317752387441497, 0.9937218650555448, 0.734298880933169]
⋮
[0.8356611821233333, 0.9685433666896768, 0.07177778586506256]
[0.8357215189457186, 0.9687394282011299, 0.07215217655460368]
[0.8360608350166773, 0.9688064835542749, 0.07224444259273748]
[0.8363380269929497, 0.9686563527743229, 0.07187609281675093]
[0.8362469897749424, 0.9684765059285643, 0.07146073353466269]
[0.8359183191543937, 0.9684329145945476, 0.07133940355770937]
[0.8356860288883594, 0.9686164623565069, 0.07169023015628033]
[0.8358384431155644, 0.9687644979956029, 0.07196671370575947]
[0.8361799095838119, 0.9687574111957825, 0.07190500575027514]
[0.8363047698281344, 0.9685635581534429, 0.07145091259050755]
[0.836100270297997, 0.968449910332632, 0.07118452254630536]
[0.836006694392986, 0.9684444145580585, 0.07116563305258432]
sol.t
1330-element Array{Float64,1}:
0.0
0.0748259153505046
0.20933786085606826
0.3735054100661737
0.5020096204826823
0.6244492914795371
0.6984427361334352
0.7504730619065657
0.7867571720102611
0.8128332086393241
0.8321016278375952
0.8471244067318621
0.8592935855891801
⋮
294.19219816611565
294.71387668708013
295.25699018406544
295.84758602400206
296.37114968486003
296.92008543869275
297.58685767529755
298.0933298151792
298.67555456940664
299.32124325626256
299.83609371027626
300.0
plot(sol.t, [[u1,u2,u3] for ((u1,u2,u3),t) in state_variables])
BoundsError: attempt to access 3-element Array{Float64,1} at index [1:1330]
Stacktrace:
[1] throw_boundserror(::Array{Float64,1}, ::Tuple{UnitRange{Int64}}) at .\abstractarray.jl:537
[2] checkbounds at .\abstractarray.jl:502 [inlined]
[3] getindex(::Array{Float64,1}, ::UnitRange{Int64}) at .\array.jl:794
[4] plotly_series_segments(::Plots.Series, ::Dict{Symbol,Any}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Tuple{Float64,Float64}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:724
[5] plotly_series(::Plots.Plot{Plots.PlotlyBackend}, ::Plots.Series) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:545
[6] plotly_series(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:862
[7] plotly_series_json at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:868 [inlined]
[8] plotly_html_body(::Plots.Plot{Plots.PlotlyBackend}, ::Nothing) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:908
[9] plotly_html_body at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:903 [inlined]
[10] html_body at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:873 [inlined]
[11] standalone_html(::Plots.Plot{Plots.PlotlyBackend}; title::String) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\web.jl:7
[12] standalone_html(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\web.jl:7
[13] _show at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\backends\plotly.jl:944 [inlined]
[14] show(::Base.GenericIOBuffer{Array{UInt8,1}}, ::MIME{Symbol("text/html")}, ::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\output.jl:215
[15] sprint(::Function, ::MIME{Symbol("text/html")}, ::Vararg{Any,N} where N; context::Nothing, sizehint::Int64) at .\strings\io.jl:105
[16] sprint at .\strings\io.jl:101 [inlined]
[17] _ijulia_display_dict(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\ijulia.jl:56
[18] display_dict(::Plots.Plot{Plots.PlotlyBackend}) at C:\Users\ibuckley\.julia\packages\Plots\hyS17\src\init.jl:73
[19] #invokelatest#1 at .\essentials.jl:712 [inlined]
[20] invokelatest at .\essentials.jl:711 [inlined]
[21] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\ibuckley\.julia\packages\IJulia\e8kqU\src\execute_request.jl:112
[22] #invokelatest#1 at .\essentials.jl:712 [inlined]
[23] invokelatest at .\essentials.jl:711 [inlined]
[24] eventloop(::ZMQ.Socket) at C:\Users\ibuckley\.julia\packages\IJulia\e8kqU\src\eventloop.jl:8
[25] (::IJulia.var"#15#18")() at .\task.jl:358
Y_output = Y0*yK(:,2)/𝜆0.*exp((𝛼+𝛽)*tK);
Error using eval Undefined function 'convert' for input arguments of type 'double'.
# Figure 3
#figure(1)
#plot(yK(:,1), yK(:,2)), hold on;
#plot(bar_omega_K, bar_lambda_K, 'ro')
#% Figure 4
#figure(2)
#yyaxis right
#plot(tK, yK(:,1),tK,yK(:,2),tK,yK(:,3))
#ylabel('\omega,\lambda,d')
#yyaxis left
#plot(tK, Y_output)
#ylabel('Y')
#title(['\omega_0 = ', num2str(omega0,txt_format), ...
# ', \lambda_0 = ', num2str(lambda0, txt_format), ...
# ', d_0 = ', num2str(d0, txt_format), ...
# ', Y_0 = ', num2str(Y0, txt_format)])
#legend('Y','\omega', '\lambda','d')
#print('keen_time_convergent','-depsc')
#omega0 = 0.75;
#lambda0 = 0.70;
#d0 = 0.1;
#Y0 = 100;
#T = 300;
#[tK,yK] = ode15s(@(t,y) keen(y), [0 T], convert([omega0, lambda0, d0]), options);
#yK = retrieve(yK);
#Y_output = Y0*yK(:,2)/lambda0.*exp((alpha+beta)*tK);
#% Figure 5
#figure(3)
#yyaxis right
#plot(tK, yK(:,1),tK,yK(:,2))
#ylabel('\omega,\lambda')
#yyaxis left
#plot(tK,yK(:,3),tK, Y_output)
#ylabel('d,Y')
#title(['\omega_0 = ', num2str(omega0,txt_format), ...
# ', \lambda_0 = ', num2str(lambda0, txt_format), ...
# ', d_0 = ', num2str(d0, txt_format), ...
# ', Y_0 = ', num2str(Y0, txt_format)])
#legend('d', 'Y','\omega','\lambda')
#print('keen_time_divergent','-depsc')
Do not run these in the notebook. MATLAB does not support functions outside of files of the same name.
### Auxiliary functions
## Instead of integrating the system in terms of omega, lambda and d, we
## decided to use:
##
## log_omega = log(omega),
## tan_lambda = tan((lambda-0.5)*pi), and
## pi_n=1-omega-r*d
## log_p = log(p)
##
## Experiments show that the numerical methods work better with
## these variables.
function new = convert(old)
# This function converts
# [omega, lambda, d, p]
# to
# [log_omega, tan_lambda, pi_n, log_p]
# It also work with the first 2 or 3 elements alone.
# each of these variables might also be a time-dependent vector
n = size(old,2); #number of variables
new = zeros(size(old));
new(:,1) = log(old(:,1));
new(:,2) = tan((old(:,2)-0.5)*pi);
if n>2
new(:,3) = 1-old(:,1)-r*old(:,3);
if n==4
new(:,4) = log(old(:,4));
end
end
end
Error: Function definition not supported in this context. Create functions in code file.
function old = retrieve(new)
# This function retrieves [omega, lambda, d, p]
# from
# [log_omega, tan_lambda, pi_n, log_p]
#
# It also work with the first 2 or 3 elements alone.
# each of these variables might also be a time-dependent vector
n = size(new,2); #number of variables
old = zeros(size(new));
old(:,1) = exp(new(:,1));
old(:,2) = atan(new(:,2))/pi+0.5;
if n>2
old(:,3) = (1-old(:,1)-new(:,3))/r;
if n==4
old(:,4) = exp(new(:,4));
end
end
end
Error: Function definition not supported in this context. Create functions in code file.
function f = keen(y)
f = zeros(3,1);
log_omega = y(1);
tan_lambda = y(2);
pi_n = y(3);
lambda = atan(tan_lambda)/pi+0.5;
omega = exp(log_omega);
g_Y = fun_kappa(pi_n)/nu-delta;
% Equation (47)
f(1) = fun_Phi(lambda)-alpha; %d(log_omega)/dt
f(2) = (1+tan_lambda^2)*pi*lambda*(g_Y-alpha-beta); %d(tan_lambda)/dt
f(3) = -omega*f(1)-r*(fun_kappa(pi_n)-pi_n)+(1-omega-pi_n)*g_Y; %d(pi_n)/dt
end
Error: Function definition not supported in this context. Create functions in code file.
tx = linspace (-8, 8, 41);
ty = tx;
[xx, yy] = meshgrid (tx, ty);
r = sqrt (xx .^ 2 + yy .^ 2) + eps;
tz = sin (r) ./ r;
mesh(tx, ty, tz)